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Abstract 

A new Particle-in-Cell (PIC) method, that conserves energy exactly, is pre- 
sented. The particle equations of motion and the Maxwell's equations are 
differenced implicitly in time by the midpoint rule and solved concurrently 
by a Jacobian-free Newton Krylov (JFNK) solver. Several tests show that 
the finite grid instability is eliminated in energy conserving PIC simulations, 
and the method correctly describes the two-stream and Weibel instabilities, 
conserving exactly the total energy. The computational time of the energy 
conserving PIC method increases linearly with the number of particles, and it 
is rather insensitive to the number of grid points and time step. The kinetic 
enslavement technique can be effectively used to reduce the problem matrix 
size and the number of JFNK solver iterations. 
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1. Introduction 

The Particle-in-Cell (PIC) method is one of the most used numerical 
methods for the solution of the collision-less kinetic equation of plasmas. The 
majority of PIC schemes has the property of conserving exactly the system 
total momentum [TJ [2] , while it does not conserve the system total energy. In 
fact typically PIC methods, that use explicit differentiation in time (explicit 
PIC schemes), tend to increase the total energy of the system by numerical 
heating [U [2] , while PIC methods, that use implicit differentiation in time 
(implicit PIC schemes), tend to decrease the total energy of the system by 
numerical cooling [3|. In the study of plasma physics instabilities, where one 
kind of energy is converted to another one, it is important to ensure that 
energy is not created spuriously by the numerical scheme in use. In fact, 
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the numerical heating introduces spurious energy that can feed erroneously 
plasma instabilities, leading to unphysical results. 

A new class of PIC methods were developed at beginning of the Seven- 
ties to address the problem of the total energy conservation. An "energy 
conserving" PIC method was first proposed by Lewis in 1970 [I]. However, 
the Lewis' "energy conserving" scheme conserves energy only in the limit 
of zero time steps, while accuracy errors preclude the possibility of having 
exact energy conservation with finite time steps [U |2] . The scheme does not 
lead to the exact energy conservation, but instead to an improved energy 
conservation at a given time step if compared to momentum conserving PIC 
methods. The "energy conserving" method was derived from the Lagrangian 
formalism, and proved to be more robust against aliasing instabilities, such 
as the finite grid instability. An in-depth analysis of the "energy conserving" 
PIC method is presented by Langdon in Ref. j5] and in textbooks jTJ 

Differently from previous "energy conserving" PIC methods, the proposed 
PIC scheme has the property of conserving exactly the total energy not only 
in the limit of zero time step, but also with finite time step. The new PIC 
method conserves the system total energy with a precision that depends 
only on the error tolerance value of the iterative solver in use. The scheme 
is not derived from the Lagrangian formalism but instead from the Newton 
approach. Energy conservation is achieved by using the midpoint integration 
rule for particle and field equations, and a convenient discretization of the 
discrete spatial operators. 

The energy conserving PIC method requires the concurrent solution of 
the equation of motion for each particle and of the field equations for the elec- 
tric and magnetic fields at each grid point, posing two challenges. First, the 
direct solution of the implicit numerical equations of PIC method implies the 
computation of a non-linear system. Non-linearity arises from the coupling 
between particles and fields variables through the interpolation functions of 
the PIC method. In the proposed PIC scheme, the non-linear equations are 
solved by a Newton Krylov solver [HI |7j. Despite the belief that the itera- 
tive solution of such equations could hardly converge [SJ, it has been proven 
that such PIC methods, that are called fully implicit, are convergent [HI [H)] . 
The second challenge is that the energy conserving PIC scheme requires the 
solution of a very large matrix whose rank is of the order of the number of 
particles (the number of particles is considerably higher than the number of 
grid points in typical PIC simulations). For this reason, implementations 
of fully implicit PIC methods are based on the matrix-free Jacobian-free 
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solvers to avoid the storage of the matrix, and Jacobian coefficients [5]. Pre- 
vious implementations of fully implicit PIC methods, as those presented in 
Refs. 0, HO] , were limited to electrostatic simulation and more importantly 
their formulation does not imply the total energy conservation. The new 
energy conserving PIC method is still based on the solution of coupled non- 
linear equations by a Jacobian-free Newton Krylov (JFNK) solver, but on the 
contrary it is formulated for the electromagnetic case and conserves exactly 
the total energy. 

This paper presents the algorithm, the properties, the implementation, 
the simulation and performance results of the energy conserving PIC method. 
It is organized as follows. Section 2 introduces the governing equations, shows 
their discretization in time and space and explains in detail the numerical 
algorithm. Sections 3,4,5 analyze the properties of the proposed method: 
the energy and momentum conservation, and the numerical stability. The 
implementation of an energy conserving PIC code is discussed in Section 
6. Section 7 presents first the results of a Maxwellian plasma simulation, 
that is robust against the finite grid instability, and then of the two-stream 
and Weibel instabilities. The computational performance of the proposed 
method, and a technique to reduce the problem matrix size are shown in 
Section 8. Finally, Section 9 concludes the paper summarizing the algorithm 
and its properties. In Appendices A and B, a skeleton version of the energy 
conserving PIC code in Matlab/Octave programming language is provided. 

2. Algorithm 

In the PIC method N s computational particles of the different n s species 
with label s mimic the real behavior of electrons and ions [11] . Each compu- 
tational particle is characterized by a position x p and a velocity v p , whose 
evolution is described by the equation of motion (here and thereafter in CGS 
units): 



electric and magnetic fields acting on the particle p and they are calculated 
by interpolation from E s and B 9 , the values of the electric and magnetic 
field on the N g grid points, through the use of the interpolation function 
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W(x g -x p ): 

E p = J2 E 9 W (*g - x p ) B p = Y^ B g W(x s - x p ). (2) 

9 9 

The Cloud-in-Cell interpolation functions [U [2] are used: 

W(x -x ) = / 1 ~ \*9-*p\/ Ax if |x 9 -x p |<Ax 
9 p ^ otherwise. 

Equations [I] are differenced in time using the implicit midpoint integration 
rule P2]: 

v™ +1 = + q s /m s At(E p + v p /c x B p ) , , 

K +1 = K + ^ P Ai ' 
where n is the time level and the bar variables are the average in time of the 
quantities, and they are defined as: 



x 



p 



x£ +1 +x£)/2 v p = (v p n+1 +v p ")/2 (5) 



E p = ]TE 9 Jy(x 9 -x p ) B p = ^B 9 iy(x 9 -x p ) (6) 

9 9 

E g = (EJ+ 1 + E£)/2 B g = (BJ+ 1 + B»)/2. (7) 

It is possible to rewrite Equations [4] in terms of v p after a series of algebraic 
manipulations [3J |9] : 

i .+ c( 4 I »B,+ g;(i,'B,)B,) 



v 



si • 



(9) 



4m,2c 2 P 

and the equation of motion becomes: 

n+l 9v — V n 

x™ +1 = x™ + v p At. 1 ; 

The evolution of the electric and magnetic fields is determined by solving the 
Maxwell's equations: 

V • E = 4vrp 

V ' B = ° (ID 

1/cdE/dt = V x B -4tt/cJ y 1 

1/cdB/dt = -V x E, 



The implicit midpoint scheme is used to discretize the Maxwell's equations. 
The Faraday's and Ampere's laws are differenced in time as follows: 



f E™ +1 -EJ = cVx B g At - 4vrJ g At = c/2V x (B™ +1 + B")At - 4ttJ s 
\ B™ +1 — = -cV x E g At = -c/2V x (E™ +1 + E")Ai. 

(12) 

The average current density J g is calculated from the particle average posi- 
tions and velocities by interpolation: 



j 9 = EE^vn x 3-^)/^> as) 

s p 

where V g is the volume of cell g. 

The discretization of spatial operators must be chosen carefully to ensure 
that the vector identities, that are valid in the continuous space, hold on the 
discrete grid also. To achieve this, the Yee's lattice [13], [H] discretization of 



the spatial operators in Equations 12 is used. Taking a uniform rectangu 



lar grid for simplicity, the different components of the electromagnetic field 
and of the current densities are calculated on the cell center (half integer 
index) and on the cell vertices (integer index) according to the Yee's lattice 
configuration: 

^9 = \^i,j+l/2,k+l/2^i+l/2,j,k^i+l/2,j,k) 

^9 = (^i+l/2J,k^iJ+l/2,h+l/2^i,,j+l/2,k+l/2) (14) 
J 9 ~~ \ u i,j+l/2,k+l/2i J i+l/2,j,k+l/2' J i+l/2,j+l/2,kJ- 

The discrete operator V is a centered difference in space (second order accu- 
rate): 

Vfjj,fc = ((fi+l/2,j,k — fi-l/2,j,k)/ Ax, (/ij+l/2,fe — fij-l/2,k)/Ay, , ^ 

(ft j,fe+l/2 - Ji,j,k-l/ 

2 )/Az). ^ 

The properties hold for the discrete operator V in the chosen spatial dis- 
cretization: 

V • V x f = V ■ (f x h) = h ■ (V x f) - f • (V x h) (16) 

The energy conserving PIC method is based on the concurrent solution 
of the coupled Equations [9] and 12 by a non-linear solver. The algorithm 
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Initialization of particle positions and velocities 
and of a self-consistent electromagnetic field 




Computational cycle 







Update of the particle new 
velocities and positions 







Calculation of the particle 
average velocities, and of the 
new values of the 
electromagnetic field by a 
non-linear solver 



J 



Figure 1: The energy conserving PIC algorithm. After the simulation has been initialized, 
the computational cycle (non-linear solver stage, and particle update) is repeated. 

is summarized in Figure [T] The simulation is initialized first, setting the 
particles positions and electromagnetic field self consistently. At each PIC 
computational cycle, Equations [9j and 12 (non-linearly coupled by Equation 



13) are solved by a JFNK solver. The governing equations have been time 
differenced with the implicit midpoint method, and solved in terms of v p , 
E n+1 , and B n+1 . The spatial operators of Equations 12 are differenced in 
space on the Yee's lattice. Once v p has been calculated with the solver, the 



new particle positions and velocities are simply updated with Equation 10 



2.1. The Electrostatic Limit 

The electrostatic formulation of the energy conserving PIC method can 
be derived by considering an unmagnetized plasma. In this case, the particle 
average velocity Equation [9] simply reduces to: 

v P = K + ^E p At = v; + ^(K +1 (*p) + E p n (x P ))At (17) 

The evolution of the electric field is determined by the Ampere's law: 



E 



n+l 



-47rJ e At. 



(18) 
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2.2. Divergence Equations 

Only the Faraday's and Ampere's equations were considered so far, while 



the two divergence equations of the system 11 were not taken in account. It 
easy to show that the equation V • B = is always satisfied if it is initially pQ . 
Moreover, the Gauss' law V • E = 4np equation is automatically satisfied by 



Equation 12, if the charge continuity equation dp/dt + V • J = holds true. 
In Particle-in-Cell methods, the charge continuity equation is not always sat- 
isfied because discrepancies between the interpolation of charge and current 
densities into the grid [lj. In fact, the definition of the current density in 



the energy conserving PIC method (Equation 13) does not satisfy the charge 
density continuity equation. In this case, the method is said not to conserve 
the charge. As pointed out in Ref. [15J, the Gauss' law can be regarded as a 
conservation principle: it is not a strictly necessary equation for describing 
the evolution of electromagnetic fields, but its violation introduces numer- 
ical errors in the simulation and might lead to unphysical behavior of the 
simulated plasma. All the energy conserving PIC simulations, based on the 
non conservative current density definition, have been first initialized solving 
the Gauss' law, ensuring there is no error due to the violation of the charge 
continuity equation initially. Then, the charge conservation has been con- 
stantly checked. The error did not grow considerably or lead to unphysical 
behavior of the simulated plasma. The numerical error can be reduced us- 
ing the pseudo current method as in Refs. [151 HE] ■ In this approach, F™ is 
defined as the violation of the Gauss' law in the cell g at the time level n, 
Fg = V • E™ — 4:7rpg , and its gradient added to the Ampere's equation: 

9 At — 9 - = cV x B g - 47rJ g + 4irdVF g , (19) 

where d is a parameter that regulates the charge conservation, and F g is 
l/2(F n+1 + F n ) (differently from Ref.pS], F is calculated at n + 1/2 and 



solved implicitly). Taking the divergence of Equation 19 



jpn+l _ rpn n+1 i n 

^--V^-(^ + V.J g ) (20) 

The left side of the equation above is the heat equation. If F g is fixed to 
zero at the boundaries initially and at each simulation time step, the error 
due to non conservation of charge diffuses away with a rate determined by 
the parameter d. Figure [2] shows the error due to the violation of Gauss' law 
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Figure 2: Maximum violation of Gauss' law divided by the maximum charge density in 
a two-stream instability simulation. The blue line represents the numerical error without 
the pseudo current correction, while the red and green lines show its evolution with a 
pseudo current correction with d = 0.1c 2 /oj pe and 0.5c 2 /oj pe respectively. 

using d equal to (no pseudo current correction), and to 0.1c 2 /u pe , 0.5c 2 /<jj pe 
in a simulation of the two-stream instability. It has been found that the 
pseudo current correction decreases the error related to the non conservation 
of charge. No major differences appeared in runs with and without the pseudo 
current correction: the plasma instabilities under study started at the same 
time and presented the same growth rate. It is clear from Figure [2] that the 
numerical error builds up slowly in the case of no pseudo current correction 
(blue line). 

3. Energy Conservation 

The discretized equations of the proposed PIC method have the property 
of conserving exactly the total energy. In fact, the variation of the magnetic 



8 



and electric fields energies (We, Wb) during the time step At is: 



1 9 

A(W E + W b ) = -J2 ((Er 1 ) 2 " (Eg) 2 + (B 9 n+1 ) 2 - (B n g ) 2 )v g (21) 



Substituting Equations 12 in the formula above: 



A(W E + W B ) = l J2g 9 (E 9 • (cV x B g - AnJ g ) - B g ■ (cV x E g At)jAtV g 

= E 3 V9 (-e 9 -J 9 -v-s 9 )a^. 

(22) 

The vector identities 16, that hold in the chosen spatial discretization, have 
been used. The first term of the equation above represents the work of the 
field on the particles, while the second term S g = c/47r(E 9 x B 9 ) is the average 
Poynting flux. Its contribution to the energy variation is zero in an isolated 



system. If the expression for the average current density 3 g (Equation 13 ) is 
substituted in the formula above, then: 

A(W E + W B ) = E^ 9 ( - E, • (E7 q s v P W(x g - 5i p )/V g ))AtV g 

= - E7 QsAtv p ■ Ef 9 E 9 W (x 9 - x p ) 

= - E7 E P S m s v p ■ (V-+ 1 - v^ 1 - E g 9 v P /c x B g ) 

= -E: s E; s V2m s ((v^) 2 -(v™) 2 ). 

(23) 

Thus the numerical scheme presented in Section 2 conserves the total energy. 
It must be pointed out that a different definition of v p in Eq. [9j using B™ 
instead of B p , still leads to the conservation energy because the magnetic 
field does no work on a charged particle. In addition, the conservation of 



energy is a consequence of the definition of the current density, Equation [13 
and PIC schemes with different techniques for the computation of current 
density to ensure charge conservation do not conserve the total energy. 



4. Momentum Conservation 

The energy conserving PIC method does not conserve momentum, as re- 
ported by Langdon in Ref. [5] for other "energy conserving" PIC schemes. The 
non conservation of momentum in energy conserving PIC schemes is due to 
spurious particle self-forces arising from the non smoothness of the current 
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Figure 3: Phase space of simulation of two cold counter-streaming electron beams with 
the energy conserving PIC method. The blue dots represent particles in the initial config- 
uration, while the red dots show the phase space at t = 13 An instability due to the 
non conservation of momentum is visible as red vortices in the phase space. 



density deposition to the grid. The relevance of the particle self-forces de- 
pends largely on the interpolation functions in use, on the number of particles 
per cell, and on the grid spacing [5]. Self-forces might trigger a macroscopic 
instability, as shown in Ref. [5]. For instance in Figure [3j two cold electron 
beams, composed by 10000 particles (in blue in the Figure [3]) are moving 
initially at opposite velocities ±0.2c in a one dimensional 2.053 c/u pe long 
periodic domain with simulation time step At = 0.5 and 64 grid cells. The 
simulation is completed with the energy conserving PIC method. An alias- 
ing instability develops phase space vortices (red dots in Figure [3]) later in 
time. A parametric study has been completed varying the grid spacing, the 
time step and v c , the characteristic velocity of the plasma (e.g. the thermal 
velocity in a Maxwellian plasma or the drift velocity in beams). It has been 
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found empirically that this aliasing instability disappears if the condition 



v c At/Ax < 1.5 (24) 

is satisfied. Thus, a modified Courant-Friederics condition, that restricts the 
maximum particle motion to one and half cell per time step, must be sat- 
isfied. The velocity v c can be connected to the fast electron scales, but no 
significant limitation arises when one desires the simulation of ion dynamics. 
It is important to note that, given a v c that is related to the electron scales 
(e.g., electron thermal velocity), the ratio At/ Ax can be still adjusted to sat- 



isfy Eq. 24 to follow the slower ion scales. In fact, such numerical condition 
does not preclude the possibility of simulations with large At, but it only 
requires that the grid spacing is chosen accordingly to satisfy the condition 



24 It must be pointed out that this constraint, that arises from the non 
conservation of momentum, is very similar to the one of the implicit moment 
PIC method [3] and for this reason the energy conserving PIC method al- 
lows simulation time steps that are as large as the one allowed by the implicit 
moment PIC method. 



5. Numerical Stability 

The numerical stability of PIC methods can be determined by studying 
the plasma numerical dispersion relation [17] , following the examples of Lang- 
don in Ref . [8] , and of Brackbill and Forslund in Ref . [3j for the electrostatic 
limit. In this approach, the particle equation of motion is linearized, and the 
electric field is assumed to have an exp(icut) dependence. The linearized equa- 
tion of motion is Fourier transformed in x, and then the perturbed charge 
density and the plasma susceptibility are calculated. Following this approach, 
the numerical dispersion of the energy conserving PIC method results: 



u pe At 2 f +co cos((w - k ■ v 



1 - (^F) 2 /o(v)^£ = 0, (25) 



sin ((oo — k • 



2 



where /o(v) is the equilibrium distribution function, u pe = \J Ami e ql / m e is 
the plasma frequency, n e is the plasma density, and k is the wave vector. In 
the case of cold plasma with /o( v ) — ^( v )> the numerical dispersion relation 
reduces to: 

tan(^)sin(^) = (^) 2 . (26) 
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The roots of the dispersion relation are always real and therefore neither 
exponential growth nor damping is present at any choice of At. For this rea- 
son, the numerical scheme is linearly unconditional stable. For comparison, 
the numerical dispersion relation of explicit PIC method for a cold plasma 
leads to exponential growth for the well known condition u pe At > 2 [1] . The 
dispersion analysis that includes grid effects is not carried out in the present 
paper, and it will part of a future work. 

In the case of implicit moment PIC methods, a 9 parameter is introduced 
in the numerical scheme to decenter in time the discretization of the field 
equations [3J. The quantities q at time level 9 are defined as (1 — 9)q n + 9q n+l . 
The numerical dispersion relation of the implicit moment PIC method has 
growing solutions for 9 < 1/2, damped solutions for 9 > 1/2, and neither 
damping nor growth for 9 = 1/2 [3j. For 9 > 0.5, the implicit moment PIC 
method damps high frequency waves, that are not resolved by the time step, 
as shown in Ref. [3j. On the contrary, in energy conserving PIC simulations 
unresolved waves are not artificially damped and hold over the simulation. 
To compare the behavior of other fully implicit PIC schemes with the energy 
conserving PIC method, a 9 parameter has been introduced in the numerical 
scheme as follows: 

t> — -\r n _l_ <?sAf -pig 

V p - V p + 2m sC y V P X + 2m s c [ yP °p)°p)l\ l + 4m2 C 2 °p > (27) 

E n+i _ e™ = cV x BjAi - An3 e g At 

gn+l _ gn = _ cV x E 6 At 

The equations above are solved concurrently by a JFNK solver. For 9 = 0.5, 
the method reduces to energy conserving PIC scheme, while for 9 > 0.5 does 
not conserve the total energy. In the latter case, the method artificially cools 
the plasma damping unresolved waves. This is clear from Figure |4j where the 
numerical dispersion relation of a plasma undergoing Weibel instability [18] 
is shown. The numerical dispersion relation has been calculated by applying 
the fast Fourier transform in space and in time to the z component of the 
magnetic field, as shown in Ref. [T^]. The Weibel instability can be seen as 
a vertical red line in both panels of Figure |4| In addition, thermal noise is 
visible as a red line, that is diagonal for small k and becomes horizontal for 
high k, only in the left panel of Figure [4] for the energy conserving PIC code 
with 9 = 0.5. The radiation field is due to aliasing errors, a feature of the 
quadratically conserving schemes [20J. On the contrary, the radiation field 
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Figure 4: Dispersion relation in the Weibel instability simulation for 9 = 0.5 (energy 
conserving PIC code) and for 9 — 0.6 (fully implicit PIC code with numerical damping). 
The Weibel instability is visible as vertical red lines in both simulations. The noise in the 
radiation field (a red line that is diagonal for small k and becomes horizontal for high k) 
is present in the energy conserving simulation (left panel, 9 = 0.5). Instead it is damped, 
and therefore not present, in the 9 = 0.6 PIC simulation (right panel). 



noise is damped in the fully implicit PIC simulation with 6 = 0.6 and not 
visible in the right panel. 

6. Implementation 

An energy conserving PIC code has been developed in the Matlab/Octave 
programming language. For implementation simplicity, the code is 1D3V pQ 
with uniform grid: the space is one dimensional while the particle velocities, 
the electric and magnetic fields have three components. The extension of the 
code to the three dimensional case is straight-forward. 

After the simulation has been initialized setting self-consistently particle 
quantities and electromagnetic fields, and ensuring that the Gauss' law is 
satisfied initially, two steps are completed at each computational cycle, as 
shown in Figure [TJ 

1. The values of the dependent variables E™ +1 , B™ +1 and v p are deter- 
mined with the JFNK solver. Given x™, v™, E™ and B" from the pre- 
vious computational cycle, v p , E™ +1 and B™ +1 are calculated solving 
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concurrently Equations [9] and 12 by a JFNK solver 
2. The new particle positions and velocities x™ +1 , v™ +1 are updated with 



Equations 10 



The JFNK method solves the non-linear system G(x) = iteratively by 
computing successive linear systems: 



<9G(x) 



dx. 



5x { = -G(xj). 



(28) 



the solution guess Xi at the iteration i of the initial non-linear system G(x) = 
is calculated as: 

x i+ i = Xi + 5x { . (29) 



The solution of the linear system 28 and the solution update (Eq. [29j com- 
pose the Newton iteration, that stops when: 



G(xj) || < e a + e r || G(x ) 



(30) 



where || ■ || is the Euclidean norm, and e a and e r are the absolute and relative 
error tolerance values. The successive linear systems 28 are solved iteratively 
by a Krylov method, the Generalized Minimal Residual (GMRes) solver [21] 
in the present study. The Krylov method convergence is adjusted at each 
Newton iteration as follow: 



dG(x) 



dx 



5x t + G(xj) || < (i II G(x 4 



(31) 



where Q is the inexact Newton parameter [21]. The number of iterations of 
the GMRes solver are called Krylov iterations. The Jacobian 9C ^^ is not 
calculated directly, but instead the Gateaux derivative is used to compute: 



dG(x) 



dx 



5xj 



lim 

e-S-0 



G(x, 



5x; 



G(xi) 



(32) 



where e is a small but finite number. Because the the Jacobian does not need 
to be formed and calculated explicitly, the method is said " Jacobian- free" . 

A guess of the particle average velocity, of the new electric and magnetic 
fields (the dependent variables) is given initially as vector x KR to the GMRes 
solver. A better estimate of these values is calculated by minimizing through 
successive solver iterations the residual r (the difference between the known 
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term b and A(^xr), the non-linear system A applied to the solution guess 
*kr) 

r = b-A(x KR ). (33) 
A function where the residual r is calculated, must be provided to the JFNK 



solver. The residual is computed by solving the Equations [9] and [12] for the 
problem unknowns, v p E™ +1 and B™ +1 , in three successive steps: 

1. Given a v p estimate in x^r, x is calculated as v p At/2, and J g is com- 



puted with Equation 13 



2. Given E™ +1 and B™ +1 estimates in x^k, E p and B p are calculated with 
Equation [6j 



3. The residual r is computed with Equations [9] and 12 

The solution of the particle equations of motion and field equations, and the 
current deposition are completed at each Krylov iteration. 

A skeleton version of the Matlab/Octave code for the electrostatic limit 
with electrons and motionless ions is presented in Appendices A and B to 
show the simplicity of the proposed PIC method. The software implementa- 
tion of the Newton Krylov GMRes solver is from the Kelley's textbook [21] 
and available at the website [22] . In all the simulations, the solver maximum 
number of Newton and Krylov iterations is set to 40, and the inexact New- 
ton parameter Q is determined by the Eisenst at- Walker formula [21]. The 
Eisenstat- Walker parameter is chosen as 0.9. In the energy conserving PIC 
method, smaller error tolerance values lead to simulations with increased 
energy conservation. This is clearly visible in Figure [5j where the energy 
history of the same simulation of the two-stream instability with different 
absolute and relative solver tolerance values is plotted. On the contrary, it 
has been found that decreasing the error tolerances does not have any effect 
in the conservation of the momentum. 



7. Simulation Results 

The energy conserving PIC codes have been tested throughly. The goal 
of these tests is first to verify the new PIC method through comparison of 
the simulation results with analytical theory, and second to show the exact 
energy conservation. In this paper, the energy conserving PIC code is first 
run in the electrostatic formulation for the problems of the finite grid and 
two-stream instabilities, and then in the fully electromagnetic case for the 
Weibel instability test. 



15 



ABS tol = 1E-7, RELtol = 1E-7 
ABS tol=5E-8,RE!_tol = 5E-8 
ABS tol = 1E-8, RELtol = 1E-8 



4.1067x10' 2 - 



4.1066x10 

4.1065x10 
LU 

"to 4.1064x10 
4.1063x10 



-2 - 



-2 L 



-2 L 



-2 _ 



4.1062x10' 2 - 




Figure 5: Comparison of energy histories in a two-stream instability simulation with 
different absolute and relative error tolerance values. Smaller tolerance errors lead to an 
increased energy conservation. 
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7.1. Finite Grid Instability 

An aliasing instability, called finite grid instability, arises in explicit mo- 
mentum conserving PIC methods when the simulation grid spacing is ap- 
proximately two times larger than Debye length (TJ EE] • The finite-grid 
instability heats non physically the plasma, until the Debye length reaches a 
value comparable to half the grid spacing. Because this instability introduces 
numerical heat in the system, it appears as a macroscopic violation of the 
energy conservation. 

To test the energy conserving electrostatic PIC method against the fi- 
nite grid instability, a Maxwellian plasma is initialized with thermal velocity 
Vthe = 0.2c in a simulation box 507ic/u pe long with 64 grid cells and 50,000 
particles. The Debye length A D = v the /u pe = 0.2c/u pe results approximately 
ten times smaller than the grid spacing Ax = 2A5c/u pe . This geometri- 
cal set-up leads to the finite grid instability if an explicit momentum con- 
serving PIC code is used. The simulation evolves over 200 computational 
cycles with time step equal to 0.5u>^}-. Therefore, the numerical constraint 
VtheAt/Ax = 0.245 < 1.5 is satisfied and numerical instability does not arise 
because of the non conservation of momentum. The absolute and relative 
solver error tolerance values are both set to 10 -8 . In addition, a simulation 
with an explicit momentum conserving PIC code, starting from the same ini- 
tial configuration, has been run to compare the results. Figure [6] represents 
the phase space (each dot represents a particle in the position-velocity space) 
of the system at t — lOOo;^ 1 for the energy conserving (red dots) and explicit 
momentum conserving (blue dots) methods. The finite grid instability in the 
simulation with the explicit PIC code is visible from the blue peaks in the 
phase space. Instead the plasma retains the initial Maxwellian distribution 
in the energy conserving PIC simulation. Figure [7] shows the total energy 
history for the two simulations. The finite grid instability produces a 2% 
energy increase in the explicit PIC simulation, while the variation with the 
energy conserving PIC method is very low, 10~ 7 %. Because the energy 
conserving PIC does not undergo finite grid instability, and does not require 
to resolve the Debye length, the proposed PIC method is well suited for sim- 
ulations with large domains and/or few grid points. A complete linear and 
non-linear analysis of the finite grid instability in the energy conserving PIC 
has not been carried out and it will be a topic of a future work. 
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Figure 6: Phase space plot of a Maxwellian plasma simulation at t = lOOo;^ 1 with 
the energy conserving (red dots) and explicit (blue dots) PIC codes. The finite grid 
instability appears as peaks of the electron distribution in the phase space in the explicit 
PIC simulation, while is not present in the energy conserving PIC simulation, where the 
plasma retains the initial Maxwellian distribution. 
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Figure 7: Total energy history of a Maxwellian plasma simulation with the energy con- 
serving PIC code (red line and left y axis) and with the explicit momentum conserving 
PIC code (blue line and right y axis). Energy is in n e m e c 2 /2 units. Because of the finite 
grid instability, the total energy in the explicit PIC simulation increases 2%. On the con- 
trary, the energy variation is limited to 10 _7 % in the case of the energy conserving PIC 
simulation. 
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7.2. Two-stream Instability 

The two-stream instability is an important phenomenon occurring in 
space physics, in the injection systems for nuclear fusion machine, and in 
particle accelerators [23] . In this problem, two electron beams move initially 
in opposite directions. The two beams extinguish as result of the beams in- 
stability A simulation of the two-stream instability has been completed with 
the energy conserving PIC code. The drift velocity of the counter-streaming 
electron cold beams is ±0.2 c; the simulation box is 2.053 c/u pe long with 
64 grid points and periodic boundaries. The number of electrons and ions 
is 200,000. The charge to mass ratio for electrons and ions is one and 1836 
respectively. The simulation time step At is 0.1 w" 1 . This set-up leads to 
v c At/Ax = 0.623 < 1.5, and thus the instability due to the non conservation 
of the total momentum does not occur. The absolute and relative tolerance 
are both set to 10~ 8 . 

The linear theory predicts a growth rate of instability for the spectral 
component k — 1 u pe /c equal to 0.35355 u pe [23] in this system configuration. 
Figure [8] shows an excellent agreement between the linear theory in red line 
and the simulation results in blue line in the linear stage of the instability. 
The energy variation for an explicit momentum conserving, and an energy 
conserving PIC code, simulating the two-stream instability are compared in 
Figure [9] The total energy in the explicit PIC code shows an approximately 
5% variation, while the variation is limited to 10~ 4 % in the energy conserving 
PIC code. An increased energy conservation can be achieved, decreasing the 
error tolerance values, as shown in Figure [5] 

7. 3. Weibel Instability 

The Weibel instability is triggered by an anisotropic temperature in the 
plasma [I~8~t 124]. The instability converts the plasma temperature anisotropy 
into magnetic field relaxing the initial particle distribution function to an 
isotropic Maxwellian. Because plasma temperature anisotropies are very 
common in laboratory, space and astrophysical plasmas, the Weibel instabil- 
ity has been thoroughly studied with PIC methods [21]. The energy conserv- 
ing PIC code has been tested in a simulation box 2irc/u pe long, with 64 grid 
points and periodic boundaries. The time step At is 0.25 u pf } and simulation 
is run for 400 computational cycles. 100,000 electrons are initialized with 
uniform distribution in space and bi-Maxwellian distribution with thermal 
velocity v t h y = 0.4, and anisotropy a = (v^ hy /v^ hx — 1) = 15. Thus, the 
condition vthxAtj Ax = 1.02 to avoid the aliasing instability due to the non 
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Figure 8: Comparison between linear theory and energy conserving PIC simulation of the 
two-stream instability. The k = \uj pe /c spectral component of the electric field (in blue 
color) grows as predicted by the linear theory (red line). The electric field is in V 1 Airn e m e c 2 
units. 
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Figure 9: Comparison of the total energy history of the explicit momentum conserving 
(blue line, right y axis) and energy conserving PIC (red line, left y axis) codes for the 
two-stream instability. Energy is in n e m e c 2 /2 units. The plot shows 5% and 1(T 4 % 
variations for the explicit and energy conserving PIC codes respectively. 
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conservation of momentum, is satisfied. The mass ratio between ions and 
electron is 1836, and ions are initialized with same temperature of electrons. 
The electric and magnetic field are initially zero. The solver absolute and 
relative error tolerance values are both set to 10~ 5 . 

The linear theory predicts a growth rate of the B z component for the spec- 
tral component k = lu pe /c equal to 0.22u pe [23J. Figure 10 compares the 
Weibel instability simulation results with the analytical calculations. Linear 
theory and simulation results are in good agreement in the linear regime of 
the instability. However, the growth of the magnetic field is not purely expo- 
nential as predicted [23], but presents an oscillation in time. This oscillation 
is caused by the radiation field noise [2H 125] • As shown in Section 5 and 
in the left panel of Figure |4j the numerical dispersion relation of the Weibel 
instability with the energy conserving PIC code shows the presence of waves 
due to the thermal noise [25]. Figure 11 shows the total energy and momen- 
tum history in a Weibel instability simulation with the energy conserving 
PIC method. The total energy variation is within 10~ 3 %. Momentum is not 
conserved and oscillates between — 0.02m P c and 0.02m P c. 



8. Performance Results 

A study of the computational performance of the proposed PIC scheme 
has been completed. A Maxwellian plasma, composed of electron and a 
background ions, is simulated with an electrostatic energy conserving PIC 
in a QAc/up e long box over a period of 1000 cycles. The other simulation 
settings vary to perform a parametric study. Table [T] shows the average of 
Newton and Krylov iterations, and the execution time for different number 
of grid points, time step, number of particles, absolute and relative error 
tolerance values. The tests have been completed on a 2.4 GHz Intel Core 
Duo, 2 GB RAM memory, Mac OS X 10.6.6 using the Matlab 7.5 and the 
code in Appendices A and B. The computational performance of the energy 
conserving weakly depends on the number of cells: increasing the number of 
grid point from 64 to 512 leads to a 20% computational time increase. The 
time step has a similar effect on performance also. The simulation with dt = 
OScUpe takes only an additional 33% computational time of the simulation 
with dt = O.lWpg 1 . Instead, the computational performance strongly depends 
on the number of particles. In fact, doubling the number of particles leads 
to doubling the computational time. In addition, increasing the number of 
particles reduces the statistical noise, and consequently the convergence in 
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Figure 10: Comparison between the spectral component k — \uj pe /c of B z and linear 
theory. Simulation and analytical results are in good agreement in the linear regime of 
the Weibcl instability. The magnetic field is in V 1 Ai:n e m e c 2 units. 
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Figure 11: Total energy and momentum histories in the Weibel instability simulation. 
Energy and momentum are in n e m e c 2 /2 and m e c units. The left red y and the right blue 
y axes correspond to the total energy and momentum respectively. The total energy is 
conserved within 10 -3 % variation, while the momentum is not conserved and oscillates 
between — 0.02m e c and 0.02m e c. 
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Table 1: Average number of Newton and Krylov (GMRes) iterations and execution time 
for different number grid points, time step, number of particles, and solver error tolerances 
(e a , e r ) for a simulation of Maxwellian plasma with an electrostatic energy conserving PIC 
code. 





dt (C^ 1 ) 


N s 


Co, 


e r 




Newton 


Krylov 


Time (s) 


64 


0.1 


100000 


io- 7 , 


10" 


-7 


3.62 


2.32 


341.86 


128 


0.1 


100000 


io- 7 , 


10" 


-7 


4.04 


2.33 


401.28 


256 


0.1 


100000 


10~ 7 , 


10" 


-7 


3.98 


2.50 


416.66 


512 


0.1 


100000 


io- 7 , 


10" 


-7 


4.02 


2.49 


413.18 


64 


0.2 


100000 


io- 7 , 


10" 


-7 


4.15 


2.45 


408.38 


64 


0.4 


100000 


io- 7 , 


10' 


-7 


4.03 


2.58 


424.98 


64 


0.8 


100000 


io- 7 , 


10" 


-7 


4.31 


2.69 


455.57 


64 


0.1 


200000 


io- 7 , 


io- 


-7 


3.25 


2.37 


617.13 


64 


0.1 


400000 


io- 7 , 


10" 


-7 


2.91 


2.39 


1,144.90 


64 


0.1 


800000 


io- 7 , 


io- 


-7 


2.79 


2.42 


2,194.00 


64 


0.1 


100000 


10 8 , 


10 


-8 


4.55 


2.34 


402.19 


64 


0.1 


100000 


10 9 , 


10 


9 


4.82 


2.42 


433.95 


64 


0.1 


100000 


IO" 10 , 


10 


10 


4.92 


2.43 


447.01 
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the Newton step. The decrease JFNK error tolerance increases the degree of 
conservation energy as shown in Figure [5] at the cost of an increased number 
of Newton iterations and computational time. 

8.1. Kinetic Enslavement 

The main disadvantage of the energy conserving PIC method is that it 
requires the solution of a very large system whose size increases with the 
number of particles. The number of particles is easily more million in a 
typical PIC simulation, leading to matrix to be inverted whose rank size 
is of the order of million. To reduce the size of this matrix, it is possible 
to use a technique, called kinetic enslavement, following Ref.[2S]- In this 
method, the unknowns of the problem are only the value of the electric 
magnetic field on the grid points at the new time level and the JFNK solver 



computes only the field Equations 12 The particle equations of motion are 



calculated by a Newton-Raphson method and embedded in the field solver 



as function evaluations [20] • Figure 12 shows the the computational cycle 



of the energy conserving PIC method with kinetic enslavement. A guess of 



the electromagnetic field (E, B in Figure 12) is provided at each Newton 
iteration. Particle positions and velocities (x, v in Figure 12) are computed 
consistently with the the electromagnetic field guess by the Newton-Raphson 
method. 

An energy conserving PIC code with kinetic enslavement technique has 



been developed to test its effectiveness. Figures [13j [14] show a comparisons 
of the solver iterations for the energy conserving PIC code with and without 
the kinetic enslavement. The plots represent the number of Newton and 
average Krylov iterations per Newton step in the simulation of the two stream 
instability with an electrostatic energy conserving PIC code with and without 
kinetic enslavement. The two stream instability is simulated for 500 cycles, 
with dt = 0.1, 64 grid points and 1000 particles. The Newton-Raphson 
convergence tolerance is set to 10 -4 , while the JFNK error tolerance values 
are set to 10~ 7 . The size of the systems is reduced from 1064 x 1064 (energy 
conserving PIC) to 64 x 64 (energy conserving PIC with kinetic enslavement). 



In addition, it is clear from Figure 13 that the use of kinetic enslavement 



technique reduces the number of Newton iterations: the average number 
of Newton iterations is 1.91 and 4.33 in the simulation with and without 
kinetic enslavement method. The number of Krylov iterations remains almost 



unchanged in the two cases (Figure 14). It must be noted that an average of 
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Figure 12: Computation cycle for the energy conserving PIC with kinetic enslavement. 
Particle positions and velocities, x, v, are calculated consistently with the fields, E, B by 
Newton-Raphson method at each JFNK iteration. 



3.95 Newton-Raphson iterations have been completed at each solver iteration. 



9. Conclusions 

The algorithm, the properties, the implementation, the simulation and 
performance results of the energy conserving Particle-in-Cell method have 
been presented. The proposed PIC method has been tested against the finite 
grid, two-stream and Weibel instabilities to prove the algorithm correctness 
and its property of conserving exactly the total energy. The method is based 
on a non conservative definition of the current density. The numerical error 
due to the violation of the Gauss' law built up slowly in all the completed 
simulations and did not affect the results. The new method does not conserve 
the momentum, and a condition on the maximum number of cells a particle 
can cross per time step, must be satisfied to avoid an aliasing instability. 
The new PIC scheme is based on the implicit discretization of the governing 
equation, and therefore results linearly unconditionally stable. 

The energy conserving PIC method is a fully implicit PIC method [TOl [9], 
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Figure 13: Comparison of number of Newton iterations in the energy conserving PIC 
code with and without kinetic enslavement for the two stream instability test. 
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Figure 14: Comparison of the average number of Krylov iterations per Newton step in 
the energy conserving PIC code with and without kinetic enslavement for the two stream 
instability test. 
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where the particle average velocities and the electromagnetic field are calcu- 
lated self-consistently in the JFNK solver to preserve the system total energy. 
The performance of the fully implicit PIC methods, and a comparison be- 
tween fully implicit and implicit moment PIC methods in terms of efficiency 
and required computational resources, have been presented in Refs.0,[HI]- It 
has been shown in this paper that the energy conserving PIC performance 
critically depends on two factors: the number of computational particles and 
the solver error tolerance values. In fact, the computational time increases 
linearly with the number of particles, and it is rather insensitive to the num- 
ber of grid points and the time step. In addition, a smaller error tolerance 
value leads to a larger number of iterations, and therefore to a larger compu- 
tational time. The use of the kinetic enslavement technique [26J reduces the 
size of problem matrix and has the beneficial effect of decreasing the number 
of Newton iterations. 
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Appendix A. ECpicES.m 



A skeleton version of the energy conserving PIC in Matlab / Octave pro- 
gramming language is here presented for the electrostatic limit (Section 2.1) 
with a plasma of electrons and motionless background ions [27]. The en- 
ergy conserving PIC code requires additional files (nsolgm.m, fdgmres.m, 
givapp.m, and dirder.m), that are available at the Kelley's textbook web- 
site [22] ■ 

After the initial parameters are set and the electric field is calculated 
solving the Gauss' law to ensure the charge conservation, the average par- 
ticle velocities and the new electric field are calculated at line 52, and the 
new particle positions and velocities are updated at lines 55 and 56 at each 
computational cycle. 

1 % Energy Conserving PIC code 

2 global L; global dx ; global NG; global DT; global N; global WP; 

3 global QM; global Q; global rho_back ; 

4 global xO ; global vO ; global EO ; 

5 

6 % simulation parameters 

r L=2*pi /3.0600; % Simulation box length 

s DT=Q.l; % time step 

9 NT=800; % number of computational cycles 

io NG=128; % number of cells 

n N=500000; % number of particles 

12 WP=1; % plasma freguency 

13 QM=— 1; % electron charge to mass ratio 

14 V0 = 0.2; % beam velocity 

is VT=0.0; % thermal velocity 

i6 tol = [IE— 7, IE— 7]; % absolute and relative error tolerance 

it dx=L/NG; % grid spacing 

is Q=WP" 2 / (QM*N/L ) ; % computational particle charge 

19 rho_back=— Q*N/L ; % background ion charge density 

20 histEnergy = []; %total energy history 

21 % two— stream instability 

22 % initial particle positions 

23 xO=linspace (0 ,L— L/N,N) '; % uniform in space 

24 % initial particle velocities 

25 vO=VT*randn(N, 1 ) ; % two counterstreaming beams 

26 pm=[l:N] pm=l-2*mod(pm, 2 ) ; vO=vO-Hpm.* VO ; 

27 

28 % Perturbation 

29 XP1=1; V1 = 0.0; mode = l; 

30 vO=vO+Vl* sin (2* pi *x0 /L*mode ) ; 
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xO=xO+XPl*(L/N)*sin(2*pi*xO/L*mode) ; 
out = (xO<0); xO(out)=xO(out)+L; 
out = (xO>=L); xO(out)=xO(out)-L; 

% calculate EO, satisfying the Gauss' Law 
% solving the Poisson equation 
p = l:N;p=[p p] ; unbones (NG— 1,1); 

Poisson=spdiags ( [un -2*un un],[-l 1] ,NG-1,NG-1); 

gl=floor(x0/dx-.5) + l; g=[gl;gl+l]; 

fr az 1 =1— abs ( xO /dx— gl + .5); fraz=[frazl;l — frazl]; 

out=(g<l);g(out)=g(out)+NG; 

out = (g>NG);g(out)=g(out)-NG; 

mat=sparse (p , g , fraz ,N,NG); 

rho=f u 1 1 ( (Q/dx) *sum(mat ) ) '+rho_back ; 

Phi=Poisson\(-rho(l:NG-l)*dx"2);Phi = [Phi ;0] ; 

EO =([Phi(NG): Phi(l:NG-l)]-[Phi(2:NG);Phi(l)])/(2*dx); 

for it=l:NT 

% start computational cycle 

xkrylov = [vO ; EO] ; 

% Newton Krylov GMRes solver 

[sol, it _hist , ierr] = nsolgm (xkrylov , 'residueEC ' ,tol ); 
v_average = sol(l:N); 

% update particle positions and velocities 
vO = 2*v_average — vO ; 
xO = xO + v_average *DT; 

% check if particle are out of the periodic boundaries 

out = (xO<0); xO(out)=xO(out)+L; 

out = (xO>=L ) ; xO ( out)=xO ( out)-L ; 

% new electric field 

EO = sol ((N+1):(N + NG)); 

% calculate the total energy 

Etot = 0.5*abs(Q)*sum(vO . " 2) + . 5 *sum(E0 . " 2 ) * dx ; 

% save the total energy 
histEnergy = [histEnergy Etot]; 
% end computational cycle 

end 



Appendix B. residueEC. m 

The Newton Krylov GMRes solver (nsolgm. m) requires the definition of 
a residue function (residueEC. m), where the discretized equations of the 
energy conserving PIC method are formulated. The following Matlab/Octave 
code presents the residue function for the energy conserving PIC code. The 
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particle average velocity equations 17 are defined at line 22, while the field 
equations 18 at line 30. 



% residual calculation for the EC PIC 
function res = residueEC ( xkrylov ) 

global L; global dx ; global NG; global DT; global N; global WP; 
global QM; global Q; global rho_back; 
global xO ; global vO ; global E0 ; 

% calculate the x at n+1/2 time level 
x_average = xO + xkrylov ( 1 :N) *DT/2 ; 

% check if particle are out of the periodic boundaries 
out = (x_average<0);x_average(out) = x_average( out)+L ; 
out = (x .average >=L ) ; x_average (out) = x_average ( out )— L ; 
% interpolation 

p=l:N;p=[p p] ; gl=floor (x_average/dx-.5) + l; g=[gl;gl + l]; 

fr az 1 =1— abs (x_average ( 1 :N) / dx— gl + .5); fraz = [( fr azl );1 — frazl ] ; 

out=(g<l);g(out)=g(out) + NG; 

out = (g>NG);g(out)=g(out)- NG; 

mat=sparse (p , g , fraz ,N,NG); 

res = zeros (N + NG, 1 ) ; 

% average velocity residual 

res ( 1:N,1) = xkrylov ( 1 :N)-vO -0.25* mat *QV1* (EO+xkrylov ((N+1):(IWG)))*DT; 
% calculate the average J 

fraz=[(frazl ) . * xkrylov ( 1 :N) ; (1 — fr az 1 ) . * xkrylov ( 1 :N) ] ; 
mat=sparse (p , g , fraz ,N,NG); 
J = full ((Q/dx)*sum(mat)) '; 

% electric field residual 

res ((N+l):(N4MJ)) = xkrylov ( (N+l ) : (N+NG)) -E0+J*DT; 
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